*-------------------------------------------------------------------------------
* TITLE: 2.figures.do
* DESCRIPTION: Code to create figures for main paper
* VERSION: MP17.0
* DATE: 08/26/21 
*-------------------------------------------------------------------------------

	/* NOTE:
		The stata command lpoly for local polynomial regression changed between 
		stata versions 15 and 16. The  optimal estimation bandwidth computed by 
		default seems to be smaller in stata 16. For consistency with earlier 
		versions of the paper and across stata versions, we set the bandwidth 
		manually such that resuls are identical to the default output of stata 15,
		whenever we use the lpoly command. The main results of the paper change 
		only marginally when the default of stata 16 or later is used instead. 
	*/

**------------------------------------------------------------------------------
** SETTINGS
**------------------------------------------------------------------------------
global Out "${Output}/figures"

* figure output specifications
local fig_w = 5
local fig_res = `fig_w' * 600 
local fig_type tif
local figure_spec as(`fig_type') width(`fig_res')
local fig_ext ".`fig_type'"

*Adjusting graphical parameters
grstyle init
grstyle set plain, horizontal 

* random seed
set seed 91621


**------------------------------------------------------------------------------
** FIGURE 1: Distribution of Productive Assets in Bangladeshi Villages: all Wealth Classes
**------------------------------------------------------------------------------


* FIGURE 1: kdensity of baseline assets ----------------------------------------

	* define sample
use "$Data/PovertyTraps_analysis.dta", clear
keep if survey == 1

	* compute axis labels 
local labels  10 100 1000 10000 100000
local xlablist 
foreach y of local labels {
	local thislabloc = ln(`y'*18.46/1000+1)
	local thislab = strofreal(`y', "%9.0fc")
	local xlablist `"`xlablist' `thislabloc' "`thislab'" "'
}

	* treatment and control pre-transfer
tw 	(kdensity Lk0 [weight = popWeight] if treat==1, lcolor (blue))  ///
	(kdensity Lk0 [weight = popWeight] if treat==0, lcolor (red) lpattern(dash)) ///
	, xtitle("Baseline productive assets" "(USD PPP, log-scale)") ///
	ytitle("density") legend(label(1 "treatment") label(2 "control")) ///
	xlabel(`xlablist') 
graph export "${Out}/kdensity_lnpAssets_baseline_pretreat_USD`fig_ext'", replace `figure_spec'

	* treatment and control post-transfer
tw 	(kdensity Lk1 [weight = popWeight] if treat==1, lcolor (blue))  ///
	(kdensity Lk1 [weight = popWeight] if treat==0, lcolor (red) lpattern(dash)) ///
	, xtitle("Baseline productive assets" "(USD PPP, log-scale)") ///
	ytitle("density") legend(label(1 "treatment") label(2 "control")) ///
	xlabel(`xlablist') 
graph export "${Out}/kdensity_lnpAssets_baseline_posttreat_USD`fig_ext'", replace `figure_spec'


* test unimodality -------------------------------------------------------------
	/*
	 * define sample
use "$Data/PovertyTraps_analysis.dta", clear
keep if survey==1 & treat==0

expand popWeight
diptest Lk0
	*/

**------------------------------------------------------------------------------
** FIGURE 2: Asset Composition and Occupation by Baseline Productive Assets
**------------------------------------------------------------------------------


* FIGURE 2 a: Composition of baseline assets -----------------------------------

	* define sample
use "$Data/PovertyTraps_analysis.dta", clear
keep if survey==1

	*generate bins of total assets
local numbins=70
su Lk0, d
local min = r(min)
local max = r(max)
local binsize = (r(max)-r(min))/`numbins'
gen bin=.
forval i=1/`numbins' {
	local thisbinLow = `min' + (`i'-1)*`binsize'
	local thisbinHigh = `thisbinLow' + `binsize'
	replace bin = `i' if Lk0>=`thisbinLow' & Lk0<`thisbinHigh'
}
bysort bin: egen freq=count(Lk0)

	* divide business assets: vehicles & tools
egen vehicles = rowtotal(asset_value6 asset_value10 asset_value12 asset_value13 )
egen tools = rowtotal(asset_value4 asset_value5 asset_value7 asset_value8 asset_value9 asset_value11)  
for var poultry_value goat_value cows_value vehicles tools land_own_total_value: gen s_X = X/(k0*1000)
collapse (mean) Lk0 (mean) s_poultry_value s_goat_value s_cows_value s_vehicles s_tools s_land_own_total_value, by(bin freq)

la var s_poultry_value "share of poultry"
la var s_goat_value "share of goats"
la var s_cows_value "share of cows"
la var s_land_own_total_value "share of land"
la var s_vehicles "share of vehicles"
la var s_tools "share of tools"

	* stacked shares
gen stack_spoultry = s_poultry_value
gen stack_sgoat = s_poultry_value + s_goat_value
gen stack_stools = s_poultry_value + s_goat_value + s_tools
gen stack_svehicles = s_poultry_value + s_goat_value + s_tools + s_vehicles
gen stack_scows = s_poultry_value + s_goat_value + s_tools + s_vehicles+ s_cows_value
egen stack_sland = rowtotal(s_*)
for var stack_*: replace X=X/stack_sland

	* compute labels
local labels  10 100 1000 10000 100000
local xlablist 
foreach y of local labels {
	local thislabloc = ln(`y'*18.46/1000+1)
	local thislab = strofreal(`y', "%9.0fc")
	local xlablist `"`xlablist' `thislabloc' "`thislab'" "'
}

	* make graph
	
		* colored version 
tw 	(area stack_sland Lk0 if Lk0<8, lwidth(0) acolor(navy) ) ///
	(area stack_scows Lk0 if Lk0<8, lwidth(0) acolor(forest_green)) ///
	(area stack_svehicles Lk0 if Lk0<8, lwidth(0) acolor(ltblue)) /// 
	(area stack_stools Lk0 if Lk0<8, lwidth(0) acolor(teal)) /// 
	(area stack_sgoat Lk0 if Lk0<8, lwidth(0)  acolor(dkorange)) ///
	(area stack_spoultry Lk0 if Lk0<8, lwidth(0) acolor(maroon)) ///
	,scheme(s2mono) graphregion(color(white)) ///
	xtitle("Baseline productive assets" "(USD PPP, log-scale)") ytitle("Share of asset type (stacked)") ///
	legend(order(6 5 4 3 2 1)) ///
	legend(label(6 "Poultry") label(5 "Goats") label(4 "Tools") label(3 "Vehicles") label(2 "Cows") label(1 "Land")) ///
	xlabel(`xlablist')
graph export "${Out}/pAssets_stacked_sharesByType_baseline_USD_color_`fig_ext'", replace `figure_spec'

		
		* grayscale version
set graph_render gdi

tw 	(area stack_sland Lk0 if Lk0<8, color(white) lwidth(.1) lcolor(black) ) ///
	(parea stack_scows Lk0 if Lk0<8, lwidth(0)  pattern(pattern8)) ///
	(parea stack_svehicles Lk0 if Lk0<8, lwidth(0)  pattern(pattern1)) /// 
	(parea stack_stools Lk0 if Lk0<8, lwidth(0)  pattern(pattern4)) /// 
	(parea stack_sgoat Lk0 if Lk0<8, lwidth(0)  pattern(pattern6)) ///
	(area stack_spoultry Lk0 if Lk0<8, lwidth(0)  acolor(gs4)) ///
	,scheme(s2mono) graphregion(color(white)) ///
	xtitle("Baseline productive assets" "(USD PPP, log-scale)") ytitle("Share of asset type (stacked)") ///
	legend(order(6 5 4 3 2 1)) ///
	legend(label(6 "Poultry") label(5 "Goats") label(4 "Tools") label(3 "Vehicles") label(2 "Cows") label(1 "Land")) ///
	xlabel(`xlablist')
graph export "${Out}/pAssets_stacked_sharesByType_baseline_USD`fig_ext'", replace `figure_spec'

set graph_render default


* FIGURE 2 b: Composition of hours worked --------------------------------------

	* define sample
use "$Data/PovertyTraps_analysis.dta", clear
keep if survey==1

	*generate bins of total assets
local numbins=70
su Lk0, d
local min = r(min)
local max = r(max)
local binsize = (r(max)-r(min))/`numbins'
macro dir
gen bin=.
forval i=1/`numbins' {
	local thisbinLow = `min' + (`i'-1)*`binsize'
	local thisbinHigh = `thisbinLow' + `binsize'
	replace bin = `i' if Lk0>=`thisbinLow' & Lk0<`thisbinHigh'
}
bysort bin: egen freq=count(Lk0)

	* drop occupations with zero hours & summarize occupations with <10 hour on avrg
drop moneylender_hours_tot rentasset_hours_tot gov_officer_hours_tot  priv_officer_hours_tot rickshaw_hours_tot 

egen hrs_others = rowtotal(hhlandcult_hours_tot factory_hours_tot foodprocess_hours_tot priv_emp_hours_tot teacher_hours_tot shopworker_hours_tot hstead_bus_hours_tot vendor_hours_tot  landrent_hours_tot shopowner_hours_tot small_bus_hours_tot priest_hours_tot politician_hours_tot gov_emp_hours_tot professional_hours_tot comhealth_hours_tot skilled_hours_tot large_bus_hours_tot fisherman_hours_tot vegetable_hours_tot)

	* define variable list 
	
		/* NOTE: 
			We exclude reported time spent with poultry  (poultry_hours_tot) 
			since raising poultry is not considered work. 
		*/
		
local occ_list hrs_others beggar_hours_tot tailor_hours_tot nonag_daylabor_hours_tot agri_daylabor_hours_tot maid_hours_tot livestock_hours_tot 
egen total_hrs = rowtotal(`occ_list')
for var `occ_list': gen s_X = X/(total_hrs)
local s_occ_list 
foreach v in `occ_list' {
	local s_occ_list `s_occ_list' s_`v' 
}

collapse (mean) Lk0 (mean) `s_occ_list', by(bin freq)
gen stack=0
foreach v in `s_occ_list' {
	replace stack = stack+`v'
	gen stack_`v' = stack
}

	* compute axis labels
local labels  10 100 1000 10000 100000
local xlablist 
foreach y of local labels {
	local thislabloc = ln(`y'*18.46/1000+1)
	local thislab = strofreal(`y', "%9.0fc")
	local xlablist `"`xlablist' `thislabloc' "`thislab'" "'
}

	* create figure
	
		* color version
tw 	(area stack_s_livestock_hours Lk0 if Lk0<8, lwidth(0.1) ) ///
	(area stack_s_maid_hours Lk0 if Lk0<8, lwidth(0) ) ///
	(area stack_s_agri_daylabor_hours Lk0 if Lk0<8, lwidth(0)) ///
	(area stack_s_nonag_daylabor_hours Lk0 if Lk0<8, lwidth(0) ) ///
	(area stack_s_tailor_hours Lk0 if Lk0<8, lwidth(0) ) ///
	(area stack_s_beggar_hours Lk0 if Lk0<8, lwidth(0) ) ///
	(area stack_s_hrs_other Lk0 if Lk0<8, lwidth(0) ) ///
	,  graphregion(color(white)) ///
	xtitle("Baseline productive assets" "(USD PPP, log-scale)") ///
	ytitle("Share of hours worked in occupation" "(stacked)") ///
	legend(label(1 "livestock") label(2 "maid") label(3 "casual agriculture") label(4 "casual non-agriculture") label(5 "tailor") label(6 "beggar") label(7 "others")) ///
	xlabel(`xlablist')
graph export "${Out}/occ_hours_share_stacked_USD_color_`fig_ext'", replace `figure_spec'
	
		* grayscale version
set graph_render gdi

tw 	(area stack_s_livestock_hours Lk0 if Lk0<8, lwidth(0.1) color(white) lcolor(black)) ///
	(parea stack_s_maid_hours Lk0 if Lk0<8, lwidth(0) pattern(pattern3)) ///
	(parea stack_s_agri_daylabor_hours Lk0 if Lk0<8, lwidth(0) pattern(pattern4)) ///
	(parea stack_s_nonag_daylabor_hours Lk0 if Lk0<8, lwidth(0) pattern(pattern5)) ///
	(parea stack_s_tailor_hours Lk0 if Lk0<8, lwidth(0) pattern(pattern6)) ///
	(parea stack_s_beggar_hours Lk0 if Lk0<8, lwidth(0) pattern(pattern8)) ///
	(parea stack_s_hrs_other Lk0 if Lk0<8, lwidth(0) pattern(pattern10)) ///
	, scheme(s2mono) graphregion(color(white)) ///
	xtitle("Baseline productive assets" "(USD PPP, log-scale)") ///
	ytitle("Share of hours worked in occupation" "(stacked)") ///
	legend(label(1 "livestock") label(2 "maid") label(3 "casual agriculture") label(4 "casual non-agriculture") label(5 "tailor") label(6 "beggar") label(7 "others")) ///
	xlabel(`xlablist')
graph export "${Out}/occ_hours_share_stacked_USD`fig_ext'", replace `figure_spec'

set graph_render default


**------------------------------------------------------------------------------
** FIGURE 4: Local Polynomial Estimates of the Transition Equation
**------------------------------------------------------------------------------

local bw_treat = 0.062
local bw_contr = 0.45
	
* FIGURE 4 a:  treatment -------------------------------------------------------

	* define sample 
use "$Data/PovertyTraps_analysis.dta", clear
keep if survey_wave==1 & stup==1 & treat==1 & Lk1<=3

	* compute axis labels 
local xlabels  400 500 600 700 800 900 1000 
local xlablist 
foreach x of local xlabels {
	local thislabloc = ln(`x'*18.46/1000+1)
	local thislab = strofreal(`x', "%9.0fc")
	local xlablist `"`xlablist' `thislabloc' "`thislab'" "'
}

local ylabels 400 600 800  1000 1200 1400  
local ylablist 
foreach y of local ylabels {
	local thislabloc = ln(`y'*18.46/1000+1)
	local thislab = strofreal(`y', "%9.0fc")
	local ylablist `"`ylablist' `thislabloc' "`thislab'" "'
}

	* create figure
tw 	(lpolyci  Lk3 Lk1, bw(`bw_treat')) ///
	(line Lk1 Lk1, lpattern(dash) lcolor(green) sort) ///
	, xtitle("Productive assets plus transfer in 2007" "(USD PPP, log-scale)") ytitle("Productive assets in 2011" "(USD PPP, log-scale)") ///
	xlabel(`xlablist')  ///
	ylabel(`ylablist', angle(horizontal)) legend(off)
graph export "${Out}/TE_lpoly_treat_USD`fig_ext'", replace `figure_spec'
	
	
	
* FIGURE 4 b: Control ----------------------------------------------------------

	* define sample
use "$Data/PovertyTraps_analysis.dta", clear
gen Lk1_placebo = Lk1 if treat==1
replace Lk1_placebo = ln(k0+(Pcows_no/1000)+1) if treat==0
keep if survey==1 & treat==0 & stup==1 & Lk1_placebo<=3

	* compute axis labels
local xlabels  10 25 50 100 250 500
local xlablist 
foreach x of local xlabels {
	local thislabloc = ln(`x'*18.46/1000+1)
	local thislab = strofreal(`x', "%9.0fc")
	local xlablist `"`xlablist' `thislabloc' "`thislab'" "'
}
local ylablist `xlablist'

	* create figure
tw 	(lpolyci Lk3 Lk1, bw(`bw_contr') )  ///
	(line Lk1 Lk1, sort lcolor(green) lpattern(dash)) ///
	, xtitle("Productive assets in 2007" "(USD PPP, log-scale)") ///
	ytitle("Productive assets in 2011" "(USD PPP, log-scale)") ///
	xlabel(`xlablist') ylabel(`ylablist') ///
	legend(off) 
		
graph export "${Out}/TE_lpoly_control_USD`fig_ext'", replace	`figure_spec'


* Find threshold from FIGURE 4 -------------------------------------------------

	* define program to find threshold 
cap program drop findthrnp
program define findthrnp, rclass
	syntax [if] [in] [, min_precision(real 0.001) bw(real 0.06)]
	tempvar grid g1 t1 before_unstableSS after_unstableSS
	
	quietly: sum Lk1 `if'
	local min_grid = r(min)
	local max_grid = r(max)
	local range = ceil((`max_grid'-`min_grid')/`min_precision')
	
	if _N < `range'{
		quietly: set obs `range'
	}
	
	quietly: gen `grid' = `min_grid' + _n * `min_precision' in 1/`range'
	
	quietly: lpoly Lk3 Lk1 `if', gen(`g1' `t1') at(`grid')  bw(`bw') nograph

	* compute threshold
	sort `g1'
	
	quietly: gen `before_unstableSS' = `g1'[_n-1] if `t1'[_n-1]-`g1'[_n-1] < 0 & `t1'-`g1' > 0 & !mi(`t1')
	quietly: gen `after_unstableSS' = `g1' if `t1'[_n-1]-`g1'[_n-1] < 0 & `t1'-`g1' > 0 & !mi(`t1')
	
	quietly count if !mi(`before_unstableSS')
	local num_remaining_SS = r(N)
	local cnt = 0
	local unstable_SSs
	
	while `num_remaining_SS' > 0 {
		local cnt = `cnt' + 1
		quietly: sum `before_unstableSS'
		local this_unstableSS_before = r(min)
		quietly: replace `before_unstableSS' = . if `before_unstableSS' == r(min)
		quietly: sum `after_unstableSS'
		local this_unstableSS_after = r(min)
		quietly: replace `after_unstableSS' = . if `after_unstableSS' == r(min)
		local this_unstableSS = round((`this_unstableSS_before'+`this_unstableSS_after')/2, `min_precision')
		local unstable_SSs `unstable_SSs' `this_unstableSS'
		quietly: count if !mi(`before_unstableSS')
		local num_remaining_SS = r(N)
	}
	
	local unstable_SS : word 1 of `unstable_SSs'
	
	return scalar thr = `unstable_SS'
end

	* define sample
use "$Data/PovertyTraps_analysis.dta", clear
keep if survey==1 & stup==1 & treat==1 & Lk1<=3

	* compute threshold	
findthrnp, min_precision(0.001) bw(`bw_treat')
disp "Threshold in log 1000BDT = " r(thr)
disp "Threshold in BDT = " (exp(r(thr))-1)*1000									
disp "Threshold in USD PPP = " = (exp(r(thr))-1)*1000 / 18.46						

	* compute standard error: 
//	bootstrap r(thr), reps(1000): findthrnp, min_precision(0.001) bw(`bw_treat')


**------------------------------------------------------------------------------
** FIGURE 5: Heterogeneity of Transition Equations 
**------------------------------------------------------------------------------

	
* define programs to find thresholds -------------------------------------------
cap program drop findthrpoly_het
program define findthrpoly_het, rclass
	syntax varlist [if] [in], by(varname)
	
	local k2: word 1 of `varlist'
	local k1: word 2 of `varlist'
	
	tempvar k1_2 k1_3
	
	gen `k1_2' = `k1'^2
	gen `k1_3' = `k1'^3
	
	reg `k2' `k1' `k1_2' `k1_3' `by' `if' `in'
	
	local a1 = _b[_cons]
	local a2 = _b[_cons] + _b[`by']
	local b = _b[`k1'] - 1
	local c = _b[`k1_2']
	local d = _b[`k1_3']
	
	mata thrpoly(`a1', `b', `c', `d')
	return scalar thr_below = r(thr)
	
	mata thrpoly(`a2', `b', `c', `d')
	return scalar thr_above = r(thr)
end

cap mata mata drop thrpoly()
mata
	void thrpoly(real scalar a, real scalar b, real scalar c, real scalar d)
	{
		ss = polyroots((a, b, c, d))
		thr = Re(ss[2])
		st_numscalar("r(thr)", thr)
	}
end

* prepare dataset --------------------------------------------------------------

	* Define earnings potential as average cow returns at village level
tempfile cowR

use "$Data/PovertyTraps_analysis.dta", clear
keep if survey==1 & cows_no>0
gen cows_no2=cows_no^2
reg livestock_inc_tot cows_no cows_no2
predict cowMReturn, res
collapse (mean) cowMReturn, by(branchid spotno)
su cowMReturn,d
gen med_cowP=r(p50)
gen aboveM_cReturn=cowMReturn>med_cowP
save `cowR', replace


	* define sample 
use "$Data/PovertyTraps_analysis.dta", clear
keep if stup==1 & treat==1 & Lk1<=3

	* merge earnings potential
merge m:1 branchid spotno using `cowR'
drop if _merge==2
drop _merge


	* prepare variables
		*  savings rate in yr 2 (2009)
gen savRate=savings/(savings+(pce_total*hhsize_adult_eq))
su savRate if survey==2 & treat & stup, d
gen aboveyr2med = savRate>r(p50) if survey==2 & savRate!=.
bysort hhid5: egen aboveM_savRate2 = max(aboveyr2med)
gen savRate2_yr = savRate if survey==2
bysort hhid5: egen savRate2 = max(savRate2_yr)
drop aboveyr2med savRate2_yr

keep if survey==1

		* BMI
replace bmi=. if bmi>50
su bmi , d
gen aboveM_bmi = bmi > r(p50) & bmi!=.

		* Anxiety
tab H_anxietyD, m
gen no_anxiety = 1-H_anxietyD

		* Any formal education
gen any_educ = resp_educ_years>0

		* non-business assets
su nonbus_assets_value, d
gen aboveM_nonbus = nonbus_assets_value > r(p50) & nonbus_assets_value!=.
	

	
* create figures --------------------------------------------------------------- 

local vars aboveM_cReturn any_educ no_anxiety aboveM_bmi aboveM_savRate2 aboveM_nonbus 
local vlabs "above median earnings potential" "any formal education" "not anxious" "above median BMI" "above median savings rate" "above median non-business assets"
 
	* compute labels
local xlabels  400 500 600 700 800 900 1000 
local xlablist 
foreach x of local xlabels {
	local thislabloc = ln(`x'*18.46/1000+1)
	local thislab = strofreal(`x', "%9.0fc")
	local xlablist `"`xlablist' `thislabloc' "`thislab'" "'
}
local ylablist `xlablist'

	* create figures
forval i=1/6 {
	
	local v : word `i' of `vars'
	local vlab : word `i' of "`vlabs'"
	local novlab = subinstr("`vlab'", "above", "below", .)
	local novlab = subinstr("`novlab'", "not ", "", .)
	local novlab = subinstr("`novlab'", "any ", "no ", .)
	
	reg Lk3 c.Lk1##c.Lk1##c.Lk1  `v' 
	cap drop pred_`v'	
	predict pred_`v'

	findthrpoly_het Lk3 Lk1  , by(`v')
	global khatLow_`v' = r(thr_above)
	global khatHigh_`v' = r(thr_below)
	
	tw 	(scatter pred_`v' Lk1  if `v'==0, connect(l i) msymbol(i O) sort lcolor(red) lpattern(longdash)	) ///
		(scatter pred_`v' Lk1  if `v'==1, connect(l i) msymbol(i O) sort lcolor(blue)) ///
		(line Lk1 Lk1, sort lpattern(shortdash) lcolor(dkgreen)) ///
		, legend(order(1 2)) legend(rows(2)) legend(label(1 "`novlab'")) ///
		legend(label(2 "`vlab'")) /// legend(label(3 "45° line")) ///
		xtitle("Productive assets post-transfer in 2007" "(USD PPP, log-scale)") ytitle("Productive assets in 2011" "(USD PPP, log-scale)") ///
		xline(${khatLow_`v'} ${khatHigh_`v'}, lpattern(dot) lcolor(black)) ///
		xlabel(`xlablist') ylabel(`ylablist', angle(horizontal))
	graph export "${Out}/Heterogeneous_TE_`v'_USD`fig_ext'", replace `figure_spec'

}

	
**------------------------------------------------------------------------------
** FIGURE 6: Difference in Differences Estimates of Long Run Dynamics in Productive
** Assets and Consumption
**------------------------------------------------------------------------------

	* define sample
use "$Data/PovertyTraps_analysis.dta", clear
keep if stup==1 & treat==1 & Lk1<=3
gen aboveT=Lk1>2.333 if Lk1!=.

	* prepare variables
gen consumption = pce_total * hhsize_adult_eq
replace pAssets=pAssets*1000

	* run regressions
local colvars
local resvars
foreach v in pAssets consumption {	
	areg  `v' i.survey_wave##i.aboveT, robust a(subdistrict)
	forval s=2/5 {
		gen `v'_b`s' = _b[`s'.survey_wave#1.aboveT]
		gen `v'_se`s' = _se[`s'.survey_wave#1.aboveT]
		local colvars `colvars' `v'_b`s' `v'_se`s'
	}
	local resvars `resvars' `v'_b `v'_se
}

	* reshape data for figures
gen one = 1
collapse `colvars', by(one)

reshape long `resvars', i(one) j(survey)
label define slab 1 "2007"  2  "2009"  3  "2011"  4  "2014"  5  "2018" 
label value survey slab

foreach v in pAssets consumption {
	gen up_`v' = `v'_b + 1.65 * `v'_se
	gen lo_`v' = `v'_b - 1.65 * `v'_se	
}

	* create figures
	
		* consumption 
twoway  (line consumption_b survey, msymbol(triangle) mcolor(red) lcolor(red) ) ///
		(rcap up_consumption lo_consumption survey , lpattern(dash) lcolor(gray) sort) ///
		,legend(off) ytitle("Total HH consumption (DiD coefficient)") /// 
		xtitle(Survey wave) xlabel(,value) 
graph export "${Out}/Longterm_coefplot_consumption`fig_ext'", replace `figure_spec'
		
		* productive assets
local j pAssets
twoway  (line pAssets_b survey, msymbol(triangle) mcolor(red) lcolor(red)) ///
		(rcap up_pAssets lo_pAssets survey , lpattern(dash) lcolor(gray) sort) ///
		, legend(off) ytitle("Productive assets (DiD coefficient)") /// 
		xtitle(Survey wave) xlabel(,value)
graph export "${Out}/Longterm_coefplot_pAssets`fig_ext'", replace `figure_spec'


**------------------------------------------------------------------------------
** FIGURE 9: Share of Ultra-Poor Households Above the Poverty Threshold as a 
** Function of the Transfer Size
**------------------------------------------------------------------------------


	* draw sample of asset shocks from control group
use "$Data/PovertyTraps_analysis.dta", clear
count if stup==1 & treat==1 & Lk1!=.
local num_Tstup = r(N)
sort hhid5
keep if treat==0 & rich==0 & deltaLk3!=.
sample `num_Tstup', count
rename deltaLk3 random_shock_control
sort random_shock_control
gen int shocks_sample_id = _n
keep shocks_sample_id random_shock_control

tempfile rshock_control
save `rshock_control'

	* define sample
use "$Data/PovertyTraps_analysis.dta", clear
keep if stup==1 & treat==1 & Lk1<3

	* randomly assign shocks from control
gen r_sort = runiform()
sort r_sort
gen shocks_sample_id = _n
drop r_sort
merge 1:1 shocks_sample_id using `rshock_control', nogen

gen Lk0_plus_rshock = Lk0 + random_shock_control
replace Lk0_plus_rshock = 0 if Lk0_plus_rshock < 0
gen k0_plus_rshock = exp(Lk0_plus_rshock)


	* compute baseline gap to threshold
gen bl_gap_rshock = exp(2.333) - k0_plus_rshock 
replace bl_gap_rshock = 0 if bl_gap_rshock<0	
	
egen mean_annual_pce = mean(pce_total)
gen bl_gap_rshock_graph = (1000*bl_gap_rshock) / mean_annual_pce


	* make figure
local obs = _N																	// extends line at top right of graph to 1.05
set obs `obs'
replace bl_gap_rshock_graph = 1.05 in `obs'

sum mean_annual_pce
local annual_pce = r(mean)
sum Pcows_no, d
local avg_trns = r(p50)
local avg_trns_share_pce = `avg_trns' / `annual_pce'
local avg_trns_share_pce_text = `avg_trns_share_pce'-0.01

distplot bl_gap_rshock_graph ///
	, ylabel(,grid) scheme(s2mono) graphregion(color(white)) title(" ") ///  
	ytitle("Share of Households" "above the poverty threshold") ///
	xtitle("Transfer value" "(share of average annual per capita consumption)") ///
	xline(`avg_trns_share_pce', lcolor(blue) lpattern(dash)) ///
	text(1.1 `avg_trns_share_pce_text' "BRAC", size(small) color(blue)) 

graph export "${Out}/share_above_t_by_transfer_nobars`fig_ext'", replace `figure_spec'
